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SUMMARY 

A computational fluid dynamics code for application to traditional incompressible flow 
problems has been developed. The method is actually a slight compressibility approach 
which takes advantage of the bulk modulus and finite sound speed of all real fluids. The 
finite element numerical analog uses a dynamic differencing scheme based, in part, on a 
variational principle for computational fluid dynamics. The code was developed in order 
to study the feasibility of damping seals for high speed turbomachinery. Preliminary seal 
analyses have been performed. 

This report was prepared by Continuum, Inc. for NASA-MSFC under contract NAS8- 
35508. 
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NOMENCLATURE 
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speed of sound 

incompressible specific heat 

total system energy 

bulk modulus of elasticity 

flux integrals, Eq. (14) 

idenity tensor 

pressure 

velocity vector 

heat transfer 

entropy 

time 

temperature 

x-, y-, z- components of velocity, respectively 

conserved variable array, Eq. (14) 

volume 

artificial compressibility 
viscosity 

kinematic viscosity 
allocation parameters, Eq. (16) 
density 
stress tensor 
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INTRODUCTION 


Leakage flow through seals in high pressure turbopumps has an important effect on the 
dynamic behavior of the pumps. The rotation of this fluid reduces the damping forces 
and generates a subsynchronous whirl instability at a speed of approximately twice the 
pump critical speed, as first shown by Black (Ref. 1). The Space Shuttle Main Engine 
(SSME) currently operates near this instability limit. Uprating of the SSME requires an 
increase in pump speeds, thus additional damping is required. Recently, von Pragenau 
(Ref. 2) suggested the use of stator roughness to produce turbulent flow in the seal. The 
flow provides additional damping for the pump, as was demonstrated analytically in Ref. 
2. From his work the impetus was provided for numerical and experimental 
investigations into the damping seals. This paper discusses the numerical modeling of 
these seals. 

A realistic analysis of the flow in the SSME turbopumps requires an unsteady liquid flow 
solver. Historically, liquid, or incompressible solutions of the Navier-Stokes equations 
have been steady state. The continuity equation is written as V • q = 0 , which 
implies an infinite sound speed, thus explicit and, consequently, unsteady methods cannot 
be obtained. This formulation also complicates the calculation of flowfield pressures, 
which are an important parameter for determining damping forces. The treatment of the 
velocity-pressure coupling then becomes of paramount importance in developing the 
numerical model for damping seals. This problem will be discussed further in the next 
section. 

The numerical procedure used in the solution method is based on the VAST (VAriational 
Solution of the Transport equations) code for compressible flows developed by Prozan 
(Refs. 3-5). The VAST code is a finite difference analysis which employs a dynamic 
differencing scheme based on the hypothesis of a variational principle for the governing 
equations. This formulation promotes the stability of the system and reduces 
computational requirements. To use the variational approach, the equation governing the 
velocity-pressure coupling must be considered when deriving the variational equation. 

Since turbulent flow is the mechanism which produces the additional damping in the seal, 
the code must include a turbulence model. There are numerous models to choose from, 
but again, when transport equation models are considered, the differencing scheme and 
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the variational equation are affected. Roughness effects must be included in any model, 
and this aspect can only be handled empirically. 

Once developed, the performance of the numerical model must be verified. Results on 
the damping seals test configuration can then be obtained and the model refined. With 
the basic model, results using LOX and LI1 2 (liquid oxygen and liquid hydrogen) can be 
compared to the analyses of Ref. 2. These cases do not make good verification problems 
because the properties, and consequently the accuracy of results, of LOX and LH 2 are 
not as well known as other fluids, such as water. The final model developed in the work 
can be used to perform parametric studies of potential SSME configurations. 
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TREATMENT OF THE VELOCITY-PRESSURE COUPLING 

The solution of the governing equations in primitive variable form for incompressible 
flows has proven to be a difficult task. The number of different equation sets used is 
nearly as large as the number of researchers working in the field. This is in sharp 
contrast to the solution of compressible gas flow in which, with the ideal gas assumption, 
the equations are well known. For incompressible and/or liquid flows, there is no 
comparable equation of state with which to close the system of equations. Density, and 
consequently pressure, are usually eliminated from the continuity equation, which 
becomes a constraint on the momentum equations. Although the pressure does appear in 
the momentum equations as spatial derivatives, neither of these equations can be solved 
for pressure without introducing directional bias (Ref. 6); nor can a transient type of 
solution be developed. Thus, from the outset one is faced with "N" equations for ”N + 1" 
unknowns. 

Early numerical solutions to incompressible flow avoided this problem by using the 
vorticity-stream function formulation, which is thoroughly treated in Ref. 6. This 
procedure suffers from many well known limitations, particularly for internal and three- 
dimensional problems (see, e.g., Refs. 6,7). With the realistic three-dimensional 
problems which are within reach by using current and future generations of 
supercomputers, the use of these equations will be severely limited. 

One of the first procedures for obtaining pressures of incompressible flows required the 
solution of a Poisson equation (Ref. 6). This equation has the form 

V 2 P = S (1) 


where S is a non-linear function of the velocities and their first three derivatives. Also, 
Neumann boundary conditions are required to solve Eq. (1). The solution to Eq. (1) must 
be obtained iteratively at each step, a very time consuming process. Another drawback 
is that the solution does not behave well when complex geometries are being studied. 
Furthermore, the solution of Eq. (1) only once at each time step does not properly 
account for the velocity-pressure coupling, as demonstrated in Ref. 8. Since it is well 
known, Eq. (1) was specifically studied for inclusion into the damping seals model; 
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however, because of the problems mentioned above, it was determined to be a non-viable 
alternative. 

An increasingly popular solution used in finite difference applications is the artificial 
compressibility developed by Chorin (Ref. 9). This formulation assumes an equation of 
state as 

P = p/6 (2) 


where 6 is the artificial compressibility. Using Eq. (2) the continuity equation becomes 
<$( — )+ V • q = 0 (3) 

'at' 


Eq. (3) allows for an explicit finite difference solution for P; however, since 6 is 
technically a relaxation parameter and has no physical basis, unsteady solutions are 
meaningless. Steady state solutions, on the other hand, have been obtained successfully 
by numerous researchers. Chang and Kwak (Ref. 10) and Kwak, et al (Ref. 11) present 
one the of the most in depth and most successful applications of the method, including 
applications to SSME problems (Ref. 12). 

A Finite Element Method (FEM) which embodies the same general principle as artificial 
compressibility is the penalty function formuation (Ref. 13). This procedure is widely 
used in pure FEM techniques, although many variants have been introduced in an attempt 
to obtain accurate results (for example, Refs. 14, 15). This further serves to illustrate 
the difficulty of obtaining pressures for incompressible flows. 

The artificial compressibility method was considered for use in developing the damping 
seals model, however, the essentially arbitrary nature of 6 in Eq. (3) creates several 
problems. Perhaps the major difficulty is in the use of the variational scheme. With no 
physical basis for 6 , the differencing procedure did not behave well. Other problems 
include the trial and error procedure required to determine 6 , and the unavailability of 
an unsteady solution. 
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There are numerous other methods for treating the velocity-pressure coupling, such as 
the transformed equations of Solomon and Szymczak (Ref. 16) or the assumed elemental 
"deviatoric stress-velocity-pressure" method of Yang and Atluri (Ref. 17). Few of these 
methods provide a set of physically accurate equations for the primitive variables. Also, 
transient problems can, at best, only be "studied" by iterating for successive steady state 
solutions, which is very costly and inefficient. Numerical damping is also necessary to 
obtain solutions in most of the available methods, including the artificial compressibility 
method. 

A recently developed procedure which attempts to overcome most of these deficiencies 
is presented in Ref. 18. This work develops the pressure solution based on a finite sound 
speed, in other words, a compressible formulation. The compressibility is used to obtain 
an equation of state, which in turn introduces the pressure into the compressible form of 
the continuity equation. The resulting equation is 

/ 9P \ - 2 

— + q • V P + pc ( V • q) = 0 (4) 

v 3 1 ' 


Eq. (4) has been successfully applied in Ref. 18 to both transient and high Reynolds 
number problems. 

The procedure to be used in developing the damping seals model is very similar to that of 
Ref. 18. The final equation differs somewhat from Eq. (4) because fewer assumptions 
will be made. The details of the analysis will be presented in the next section. 
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THEORY 


Basic Equations 


In the development of the compressible VAST code, the integral form of the equations of 
motion are used. Several advantages of this formulation occur, as pointed out in Ref. 5. 
The fundamental equations are: 


continuity: 


/// — dV + // pq • ds = 0 

cv 3t cs 


(5) 


momentum: 


III 


9 Pq 

dV + 

St 


// C ( Pq ) q + ( I P- x ) ] *ds 

cs 


0 


( 6 ) 


energy: 


/// 

cv 


SpE 

dV + 

St 


// [pEq + (IP-T) 

CS 


q+Q]*ds = 0 


(7) 


For this analysis, as with all incompressible methods, the assumption is made that P 4 P 
(T), hence the energy equation (Eq. 7) decouples from the other equations. The equation 
is listed here since, if heat transfer or thermal data were desired, Eq. (7) could be solved 
independently for temperature. This equation will not be used in the remainder of the 
analysis as only isothermal problems will be considered. Since liquid flows are to be 
studied, it should also be pointed out that gravity, if assumed to be important, is easily 
included in Eq. (6) as a source term. 
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Pressure Solution 


As previously discussed, the current approach to the pressure solution is to take 
advantage of the finite sound speed in all fluids, liquid or gas. The procedure is more 
correctly called a slight compressibility method, rather than an incompressible method. 

An equation of state is assumed to exist in the form 

P = P(p) (8) 


The compressiblity of any fluid is expressed by the bulk modulus of elasticity, 


dP 



P 


(9) 


and the finite acoustic speed in a liquid is given by 

a 2 = = — 

P dp 


The equation of state may now be written as 
dP = a 2 dp 


( 10 ) 


(ID 


To derive the pressure equation, Eq. (11) is primarily used in differential form, however, 
a reference state must also be assumed for integration purposes. This reference state, 
P r and p^ , must be assumed to be the usual theoretical incompressible values of P and 
p for the fluid being considered. Substitution of Eq. (11) into Eq. (5) now yields the 
pressure equation. 
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9P 

— dV + 

at 


j j Pq *ds + (p a 2 
cs r 



q *ds = 0 


( 12 ) 


The convective part of Eq. (12) is written as two terms to facilitate discussion of some 
salient features of the method. First, it should be noted that as the sound speed, a , 
approaches infinity (the true incompressible case) Eq. (12) reduces to the standard 
incompressible form of continuity. Another feature is that the artificial compressibility 
method can be viewed as a limiting case of Eq. (12). If the contribution of the pressure 
to the convective terms is discarded and Pj.a 2 = 1/5, then Eq. (12) reduces to the 
integral form of Eq. (3). Perhaps the most important feature is the physical reasoning 
used to obtain Eq. (12). In addition to more correctly describing the physics of the 
flowfield, this method allows the use of a variationally based dynamic differencing 
scheme. 

Eq. (11) could, and perhaps should also be included in the momentum equations. The 
error introduced by using p = in Eq. (6) is essentially zero for steady state 
solutions, which can be viewed as incompressible solutions. For unsteady solutions, the 
error is of order Ap/pa^ vs. p and, since a is large, this error is assumed to be 
negligible. 


Numerical Analog 

The numerical integration used to solve Eqs. (6) and (12) is based in part on the 
variational procedure developed by Prozan in Refs. 3-5. This procedure will be 
summarized below, and the application to the current work will then be discussed. 


The compressible VAST transport equations may be generalized as 


9 

9 1 



U. 


dV 


n — 1 , . . . , 5 


where 


U 



F = 


// 


pq • ds_ _ 

[Pqq + (IP - t) ] • _ds 

[pEq + (TP-t) • q + Q ] • ds 


(13) 


(14) 
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Now, the domain of integration is subdivided into E finite elements and Eqs. (13) become 


E K 

l ( AU AV/At ) = - Z F 


(15) 


A functional distribution of the variables is assumed over the surfaces of the elements 
and the flux terms F n>e are evaluated using values from the previous time step. 


The ( AU n AV/At) g represent the total accumulation of the nth conserved quantity 

within element e. To obtain the conserved variables at a new time step, nodes are placed 

at the corners of each element and the ( AU AV/At) are allocated to the various 

n e 

nodes. Thus, the conserved quantities at the nodes are determined by assembling the 
contributions from its surrounding elements. This accumulation can be written 


(AU n AV/At). 


k 

E £ • F 

0 _i n , e , j n, e 


(16) 


where j is the node number, k is the total number of elements surrounding node j, and the 

£ . are the allocation parameters, 

n , e , j 

The assembled equation, Eq. (16) is analogous to a finite difference expression in which 
the spatial transformations are numerically embedded in the analog. It may be 
interpreted as a general form of the finite difference scheme since different finite 
difference algorithms can be derived with selected allocation parameters. The VAST 
differencing procedure differs from other schemes in that it does not dictate a fixed 

5 . throughout the course of integration, but rather, changes the parameters 

dynamically in both time and space according to the variational principle. This principle 
requires that the rate of entropy production be maximized. The VAST differencing 
scheme thus uses the transport equations as equality constraints and dynamically 
determines the allocation parameters to achieve maximum stability of the system. 

To apply the VAST numerical analog to the current investigation, a procedure for 
determining the allocation parameters similar to the procedure of the Ref. 5 must be 
developed. A first attempt of this development was presented in Ref. 19. This method 
was not entirely satisfactory, and further research was performed. The entropy 
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functional of Ref. 19 was further manipulated and the following derivatives were 
obtained: 


9 P S r 1 

= I —rr c 1 nT 

9P L a 


1 

c + - 
2 



(17) 


where c is the specific heat, 
9 ps u 

9 pu T 

9 ps v 

9 pv T 

9 ps w 

9pw T 


(18a) 


(18b) 


(18c) 


Using Eqs. (17) and (18), the allocation parameters of Eq. (16) are determined in the same 
manner presented in Ref. 5. This procedure provided good stability for many cases, 
however, some problems were encountered with the pressure equation. Eventually, a 
centered difference scheme was used for the pressure equation for all time and all nodes 
while continuing to use Eqs. (18) for the momentum equations. This method has, thus far, 
produced satisfactory results on all the problems attempted, although some damping of 
Eq. (12) is necessary. 


Turbulence Model 


To accurately model the damping seals problem, a turbulence model must be included in 
the code. Turbulent flow is an important parameter in the performance of the damping 
seals, and roughness is used to produce this flow. Hence, the model must ultimately 
include roughness effects. All turbulence models in current use have a considerable 
amount of empiricism included. Even the more advanced Reynolds stress transport 
equations require several empirical constants for closure. Roughness effects in any 
model are, at best, accounted for by these empirical constants, thus it was decided to use 
an algebraic model in the code development. 
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To develop a turbulence model for use in the VAST codes, special wall element functions 
were developed to avoid using an excessive number of grid points near the wall. This 
model was discussed in Ref. 19, but will be included here for completeness. For smooth 
walls and fully developed pipe flow, the following empirical logarithmic profiles are 
valid. 


+ + 
u = y 


+ 

y < 5 

(19a) 

+ 

. + . 

+ 


u = -3.05 

+ 5.0 1 n (y ) , 

5 < y <30 

(19b) 

+ 

+ V 

+ 


u =5.5 + 

2.5 1 n (y ) , 

y > 30 

(19c) 


where 

+ * 

u = u / u 

u* = Vt / p 

w 

+ * . 
y = u y / v 

and y is the distance from the wall. From Eq. (19a), the velocity gradient at the wall is 



( u* ) ^ / v 


(29) 


In terms of real distance from the wall, Eq. (19c) represents most of the boundary layer, 
therefore the following computational procedure is employed. A fictitous wall is 
assumed to be 0.0005 feet (approximately y + = 30) away from the real wall and it is 
assumed that no flow occurs between these two walls (m=0) . If a slip boundary 

condition is used, u is determined from the momentum equations, and u* at 0.0005 feet is 
determined from Eq. (19c). Since this equation is not explicit in u*, an empirical 
approximation is made: 

u* = 0. 1662529 u 0 - 867325 / ( y / v) 0 * 1 3 2 6 7 5 ( 2 1 
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Eqs. (20) and (21) now give the velocity gradient and shear stress at the wall. An eddy 
viscosity is used to determine both the local shear stress and the variation of this stress 
with distance from the wall, y. 

P T = 0.07 u prt + p (22) 


where r is the radius for pipe flows, and 

£ = (y/. 3r ) , 0 < (y/r ) < . 3 
£ = 1.0 , (y/r) > .3 


To include roughness in the turbulence model outlined above, the empircal constants of 
Eqs. (19) - (22) must be changed accordingly. Ref. 20 presents the information necessary 
for these changes. 
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RESULTS 

In order to verify the performance of the model during development, several simple two- 
dimensional test cases were run. Results from these problems were presented in Ref. 
19. Those results were obtained on the Continuum CM-1000 workstation (Intel 8086 
microprocessor). One important conclusion that came from the early runs was the long 
run times associated with the CFL constraint. These run times coupled with the 
geometric necessity to use a large number of grid points to model damping seals, 
dictated a shift in the development effort to a mainframe computer. 

Vectorization 


The computer code was installed on a CRAY-1S vector processor and extended to handle 
three-dimensional geometries. In order to take advantage of the available computing 
power, this version of the code was vectorized. The vector FORTRAN used is identical 
to the language to be used on the MSFC CRAY X-MP. This will allow for installation and 
efficient operation of the code at NASA MSFC. Once the mainframe programming was 
completed the two-dimensional problems presented in Ref. 18 were repeated. Results 
from the mainframe code were in agreement with the previous results, thus other 
problems with similarities to the damping seals geometry were attempted in order to 
refine the model. 

Rotating Cylinders 

Several cases involving two concentric rotating cylinders were studied. This two- 
dimensional flow is the Couette flow described in Ref. 2 for the damping seals. Both the 
differencing scheme and boundary conditions for use in the liquid code were examined 
using this configuration. This study led to the decision to use Eq. (18) with a centered 
scheme for pressure in the differencing algorithm. The boundary conditions used on the 
moving wall were a no-slip condition with a table lookup algorithm to maintain a 
constant rotational velocity at the wall. The boundary condition must be modified to 
allow a slip wall when the turbulence model is used, however, the two-dimensional results 
were obtained for a laminar viscosity. 

Velocity vectors for a 10° slice of the cylinders are shown in Figs. 1 and 2 for both air 
and water as fluids. The surface velocity of the inner cylinder is 200 ft/sec while the 
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outer cylinder is stationary. These results are compared to the exact solution in Fig. 3. 
The agreement is seen to be quite good. The Reynolds numbers of these problems were 
of order 10 ® (air) and 10® (water). This range is two to three orders of magnitude higher 
than any problems attempted by the majority of existing incompressible codes. One 
reason these large Reynolds numbers can be studied is that the integral formulation, Eqs. 
(5-7), (12), eliminates cell Reynolds number problems. 

The cylinder problem was also run using the compressible VAST code for comparison 
purposes. These results are also shown in Fig. 3. They do not agree as well with the 
exact solution. This disagreement can be traced to the behavior of the pressure 
solution. The initial condition for this problem assumed a constant pressure of 2116 
psf. The pressure after 44,000 iterations is shown in Fig. 4. (For the exact solution, the 
outer cylinder pressure is 2119.33 psf). The constant initial condition caused waves to be 
generated in the flowfield. These waves were observed to "bounce” back and forth 
between the cylinders. The no-slip walls are non-porous to these waves, and the 
absolutely conservative nature of the numerical analog (Ref. 5) coupled with the absence 
of explicit artificial viscosity do not allow the waves to damp out in this closed system 
other than through real viscous dissipation. This was verified by running the same 
problem with the exact solution as an initial condition. No waves were generated, and 
the pressure essentially did not vary from the exact solution. 

Similar pressure waves can be observed in the incompressible code. This is clearly 
illustrated by the two waves seen in Fig. 5 for the cylinder problem with water. This 
result together with the velocities (Fig. 3) tends to support the contention that the 
pressure does not greatly affect the velocity solution for essentially incompressible 
flows. This idea is important to solution procedures which eliminate the pressure 
entirely (for example the vorticity-stream function formulation) or decouple the pressure 
from the momentum equations. This small dependence is also important to the 
performance of the artificial compressibility method, as the pressures computed during 
the relaxation are entirely unrealistic. In contrast, Fig. 3 illustrates that compressible 
solutions, even at very low speeds, are sensitive to pressure effects. 

In an attempt to reduce the propagation of pressure waves in a closed system, wall 
boundary conditions were examined which would allow these waves to pass through, such 
as occurs at inlets for an open system. This effort was not entirely successful, thus some 
damping is used in the pressure differencing algorithm. 
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Fig. 4. Pressures for Cylinders with Air from Compressible VAST Code. 
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Exact Solution 


O Step 50,000, AT= .0040 sec 



Pressure, psf 


Fig. 5. Radial Pressure Profile for Rotating Cylinders 
with Water. 
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The cylinder results shown in Figs. 1-5 were obtained using a mesh of 15 x 241 (radial x 
circumferential) nodes. Thus, an element side is 1.5° of arc in length. To determine the 
effect of mesh size, the same configuration (with air as the fluid) was repeated using a 
15 x 91 mesh, or 4° length of each element. This solution showed poor agreement with 
the exact solution, proving that good circumferential resolution is required at least when 
the flow is circumferential in nature. 

Three-Dimensional Results 


After refining the model based on the two-dimensional rotating cylinder studies, a three- 
dimensional damping seals configuration was modeled and preliminary results obtained. 
The grid for the seal, based on the test configuration to be run, is shown in Fig. 6. The 
gap between the cylinders in magnified in this drawing; actual dimensions are 0.02 inch 
gap, 1.8 inch rotor diameter and 1.8 inch axial length. The grid resolution is 11x105x26, 
radial, circumferential and axial, respectively, for a total of 30,030 nodes. This 
resolution produces typical element length ratios illustrated in Fig. 6. This element 
shape in conjunction with the two-dimensional mesh study, immediately suggest that 
more nodes should be used. However, available computer storage on the CRAY-1S will 
allow 31,000 nodes maximum for this problem. Expanded inlets and outlets to the seal 
should also be included in an in-depth study, but again, this requires additional nodes. 
This node limitation will be alleviated by use of the MSFC CRAY X-MP, which will have 
more than twice the storage available on the CRAY-1S used in this study. It should be 
pointed out that the large number of nodes required is as much of a physical geometry 
requirement, as it is a numerical requirement. 

The fluid in the three-dimensional study is water ( p r =l. 89936 slugs/ft^, P r =2116 psf) at 
614 °R, and the turbulence model was used. The upstream axial boundary condition was a 
fixed inlet with the velocity profiles (both axial and radial) based on a 1/7 power law 
profile. The downstream axial boundary condition is a constant integrated mass flow in 
the axial direction. The rotor speed is 345.6 ft/sec and the axial static pressure drop 
through the seal is 2000 psi. 

Preliminary results for this problem are shown in Figs. 7-9. This is at an early iteration 
step, not at a steady state. Further results have not been obtained due to needed 
problem refinement that can only be provided by the CRAY X-MP. Furthermore, 
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experimental data is not yet available for comparisons. Fig. 7 shows the velocity vectors 
on the rotor surface for a 105° slice which has been "unrolled”. The initial condition was 
the 345.6 ft/sec speed circumferentially, with 0.0 ft/sec axial flow at the wall. The 
development of the axial flow based on the "slip condition" in the turbulence model is 
seen to be quite good. Since the upstream inlet (top of figure) is held at zero axial 
velocity although there is a circumferential component (which does not show on the scale 
used), the flow must develop along the length of the rotor. This development is also 
proceeding well. This result suggests that, to get the proper characteristics in the seal, 
the upstream boundary condition should be changed, or else some type of inlet extension 
should be included, preferably the latter choice to more correctly model the test 
configuration. 

Studying the circumferential direction at a fixed axial location in Fig. 7, some 
unexpected variations in the velocity vectors are seen. An examination of the full 360° 
data shows that this pattern repeats at exactly 90° intervals. This suggests a dependence 
on the grid resolution and the element shape, as discussed on the two-dimensional 
problem, and seen in Fig. 6. 

At an axial station in the fully developed flow region of Fig. 7 (at approximately 75 % of 
the seal length) velocities and pressures are plotted in Figs. 8 and 9, respectively. The 
axial component of the velocity is not illustrated in Fig. 8. Once again, the velocities at 
the walls are seen to be different from the cylinder velocities. This is expected due to 
the turbulence model, since this is actually the velocity at y + = 30 away from the walls. 
The turbulent profile is seen here to be developing quite well. Figure 9 shows the 
pressures at the same surfaces. The pressures are behaving well, with no large waves 
propagating at this location. 

Computational Requirements 

From the two- and three-dimensional results just discussed, mainframe computer 
requirements can be established for the code. As previously mentioned, storage 
requirements are large and indeed, the limit of 31,000 nodes for three-dimensional, 
turbulent problems is too small. Considerably larger numbers of nodes can be used for 
three-dimensional laminar, or any two-dimensional problems. This nodal limit will more 
than double on a CRAY X-MP. CPU times for a CRAY-1S computer have been 
established as follows: 
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two-dimensional, viscous, laminar problems, 1.8 x 10 -5 sec/node/step 
three-dimensional, viscous, turbulent problems, 8.9 x 10" 5 sec/node/step 

These times will be approximately a factor of 10 faster on a CRAY X-MP. Comparison 
of these resource requirements with other CFD codes show that they are very similar, 
and perhaps even smaller than existing methods (see for example Ref. 12). 
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CONCLUSIONS 


One of the principle objectives of the research work reported herein was to produce a 
high Reynolds number, explicit CFD code for liquid and/or incompressible flow regimes. 
Considerable effort was expended on this task with good results. The code which has 
been developed has been applied to several problems, with results comparable to any 
existing code, including the much studied artificial compressibility codes. Computer 
requirements necessary to generate solutions is also comparable to or less than 
requirements of other codes. However, flow problems at much higher Reynolds numbers 
than exist in the majority of the literature have been presented. There is much room for 
improvement in the model to further increase its utility. In particular, the proper 
accounting for pressure in the variational approach, and consequently the differencing 
scheme, would enhance the stability and overall performance of the method. 

Two-dimensional problems similar in character to the damping seals have been studied, 
and the computer model refined based on these results. Preliminary studies of a three- 
dimensional seal model have also been made. However, several comparisons with other 
seal analyses have not been made. In particular, the anticipated experimental data which 
is important to the model verification has not become available during the course of this 
study. Liquid oxygen and liquid hydrogen studies(which are presented in the bulk flow 
analysis of Ref. 2) are not good problems for use in model refinement due to the 
extreme physical characteristics of the fluids. Such refinement should be made prior to 
performing parametric studies because of the large amount of computer resources 
necessary to perform these studies. 

To improve the efficiency of the computer code, so that the experimental and 
parametric studies can be performed once the data becomes available, the code has been 
vectorized and developed to run on the CRAY X-MP to be installed at MSFC. This will 
also help provide NASA in-house capabilities to model turbomachinery. 
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